#include <src/pcps.h>
#include <src/input.h>
#include <src/operator.h>
#include <src/hamiltonian.h>
#include <src/wavefunction.h>
#include <src/correlator.h>
#include <src/cluster.h>
#include <src/expectation_value.h>
#include <src/blas_lapack.h>
#include <src/parallel.h>

//-------------------------------------------------------------------------------
// pcps::projected_nr_update -- updates the correlators via a Newton-Raphson step
//-------------------------------------------------------------------------------

template <class S>
void pcps::projected_nr_update(const pcps::Input & userinp,
                               const pcps::Hamiltonian<S> & ham,
                               pcps::Wavefunction<S> & wfn,
                               const int perturbation_level,
                               const double damping) {

  // not yet implemented for complex numbers
  if (typeid(S) != typeid(double))
    throw pcps::Exception("pcps::projected_nr_update is only implemented for real numbers.");

  // define the worker class that will compute the residual and jacobian precursers in parallel
  class Local_Worker : public pcps::Parallel_Worker_Base<S> {

    public:

      void do_work(const pcps::Input & userinp,
                   const pcps::Hamiltonian<S> & ham,
                   const pcps::Wavefunction<S> & wfn,
                   const int perturbation_level,
                   const boost::shared_array<S> & output_array,
                   const int array_size,
                   const int my_min,
                   const int my_max,
                   const std::vector<int> & slaves,
                   int & slaves_used,
                   int & jobs_in_progress)
      {

        // define function object to print progress of jacobian construction
        struct PrintProgress{
          void operator()(const int total, const int finished, int & recorded, const pcps::Stopwatch & sw) {
            const double old_completion = double(recorded) / total;
            const double new_completion = double(finished) / total;
            if ( int(100*new_completion) > int(100*old_completion) ) {
              std::cout << boost::format("%4i%%  %10.4f elapsed seconds") % int(100 * new_completion) % sw.elapsed_seconds() << std::endl;
              recorded = finished;
            }
          }
        };
        PrintProgress print_progress;

        // get MPI info
        const MPI::Comm & comm = MPI::COMM_WORLD;
        const int nproc = comm.Get_size();
        const int myrank = comm.Get_rank();

        // determine the number of variables in the wavefunction
        const int nvars = wfn.nvars();

        // The output array's size is nvars*(nvars+3).
        // The first nvars*nvars elements are the jacobian precursers
        // The next nvars elements are the energy derivative precursers
        // The next nvars elements are the residual precursers
        // The last nvars elements are the overlap matrix elements

        // determine how many expectation values on this processor are required to build the jacobian
        int finished_ev = 0;
        int recorded_ev = 0;
        int total_ev = nvars * ( my_max - my_min + 1 );
        for (typename std::vector<pcps::Correlator<S> >::const_iterator corr = wfn.correlators().begin(); corr != wfn.correlators().end(); corr++)
          if (!corr->fixed())
            total_ev += 2 * corr->nelements() * ( my_max - my_min ) * ( 1 + nvars );

        // prepare a copy of the identity operator
        const typename pcps::OpBase<S>::const_ptr identity_op(new pcps::OpIdentity<S>(pcps::unity<S>()));

        // prepare an empty set of correlators to indicate that none should be excluded during expectation value calcualtions
        typename pcps::Correlator<S>::Set excluded;

        pcps::Stopwatch sw; sw.start();
        double completion_counter = 0.0;

        // loop over the correlators that form the residual
        for (int r_ind = 0; r_ind < wfn.correlators().size(); r_ind++) {

          // get a reference to the residual correlator
          const pcps::Correlator<S> & res_corr = wfn.correlators().at(r_ind);

          // fixed correlators are not included
          if (res_corr.fixed()) continue;

          // loop over the correlator's occupations
          for (pcps::Cluster res_occs(userinp.nsites(), userinp.configs_per_site(), res_corr.site_vec()); !res_occs.finished(); res_occs++) {

            // get the residual element's compound index
            const int cmpd_res_index = wfn.corr_offsets()[r_ind] + res_occs.cmpd();

            // compute the contribution to the energy derivative element
            const int energy_der_pos = nvars*nvars + cmpd_res_index;

            // numerator (bra_occs = none, ket_occs = res_occs)
            S coeff = pcps::unity<S>() / res_corr[res_occs.cmpd()];
            for (int i = my_min; i < my_max; i++) {
              pcps::expectation_value<S>(slaves, slaves_used, jobs_in_progress, output_array.get(), energy_der_pos, coeff,
                                              userinp, wfn, ham.ops().at(i), perturbation_level, pcps::Cluster(), res_occs, excluded);
              if (myrank == 0 && userinp.verbose_print()) print_progress(total_ev, ++finished_ev, recorded_ev, sw);
            }

            // denominator (bra_occs = res_occs, ket_occs = none)
            coeff = -pcps::unity<S>() / res_corr[res_occs.cmpd()];
            for (int i = my_min; i < my_max; i++) {
              pcps::expectation_value<S>(slaves, slaves_used, jobs_in_progress, output_array.get(), energy_der_pos, coeff,
                                              userinp, wfn, ham.ops().at(i), perturbation_level, res_occs, pcps::Cluster(), excluded);
              if (myrank == 0 && userinp.verbose_print()) print_progress(total_ev, ++finished_ev, recorded_ev, sw);
            }

            // residual, overlap, and jacobian precursers are not computed for images of translationally invariant correlators
            if (res_corr.ti()) continue;

            // compute the residual precurser element
            const int res_precurser_pos = nvars*nvars + nvars + cmpd_res_index;
            for (int i = my_min; i < my_max; i++) {
              pcps::expectation_value<S>(slaves, slaves_used, jobs_in_progress, output_array.get(), res_precurser_pos,
                                              pcps::unity<S>(), userinp, wfn, ham.ops().at(i), perturbation_level,
                                              res_occs, pcps::Cluster(), excluded);
              if (myrank == 0 && userinp.verbose_print()) print_progress(total_ev, ++finished_ev, recorded_ev, sw);
            }

            // compute the overlap matrix element
            if (myrank == 0) {
              const int overlap_pos = nvars*nvars + 2*nvars + cmpd_res_index;
              pcps::expectation_value<S>(slaves, slaves_used, jobs_in_progress, output_array.get(), overlap_pos,
                                              pcps::unity<S>(), userinp, wfn, identity_op, perturbation_level,
                                              res_occs, pcps::Cluster(), excluded);
              if (myrank == 0 && userinp.verbose_print()) print_progress(total_ev, ++finished_ev, recorded_ev, sw);
            }

            // loop over the correlators that we differentiate with respect to
            for (int d_ind = 0; d_ind < wfn.correlators().size(); d_ind++) {

              // get a reference to the derivative correlator
              const pcps::Correlator<S> & der_corr = wfn.correlators().at(d_ind);

              // fixed correlators are not included
              if (der_corr.fixed()) continue;

              // loop over the derivative correlator's occupations
              for (pcps::Cluster der_occs(userinp.nsites(), userinp.configs_per_site(), der_corr.site_vec()); !der_occs.finished(); der_occs++) {
                
                // get the derivative element's compound index
                const int cmpd_der_index = wfn.corr_offsets()[d_ind] + der_occs.cmpd();

                // get the jacobian element's compound index
                const int cmpd_jac_index = cmpd_res_index * wfn.nvars() + cmpd_der_index;

                // compute the contribution from each hamiltonian operator
                for (int i = my_min; i < my_max; i++) {

                  // check if this operator contributes
                  bool contributes = false;
                  for (int j = 0; j < ham.ops().at(i)->sites().size(); j++)
                    if (    ham.ops().at(i)->offsets().at(j) != 0
                         && std::count(der_corr.site_vec().begin(), der_corr.site_vec().end(), ham.ops().at(i)->sites().at(j)) > 0)
                      contributes = true;

                  // if the operator does not contribute, skip it
                  if (!contributes) {
                    ++finished_ev;
                    if (myrank == 0 && userinp.verbose_print()) print_progress(total_ev, ++finished_ev, recorded_ev, sw);
                    continue;
                  }

                  // compute the numerator contribution (bra_occs = res_occs, ket_occs = der_occs)
                  coeff = pcps::unity<S>() / der_corr[der_occs.cmpd()];
                  pcps::expectation_value<S>(slaves, slaves_used, jobs_in_progress, output_array.get(), cmpd_jac_index,
                                                  coeff, userinp, wfn, ham.ops().at(i), perturbation_level,
                                                  res_occs, der_occs, excluded);
                  if (myrank == 0 && userinp.verbose_print()) print_progress(total_ev, ++finished_ev, recorded_ev, sw);

                  // compute the denominator contribution (bra_occs = res_occs + der_occs, ket_occs = none)
                  std::set<int> bra_occ_sites;
                  bra_occ_sites.insert(res_occs.sites().begin(), res_occs.sites().end());
                  bra_occ_sites.insert(der_occs.sites().begin(), der_occs.sites().end());
                  pcps::Cluster bra_occs(userinp.nsites(), userinp.configs_per_site(), bra_occ_sites);
                  bool compatible = true;
                  for (std::set<int>::const_iterator site = bra_occ_sites.begin(); site != bra_occ_sites.end() && compatible; site++) {
                    if (    std::count(res_corr.site_vec().begin(), res_corr.site_vec().end(), *site) > 0
                         && std::count(der_corr.site_vec().begin(), der_corr.site_vec().end(), *site) > 0 )
                      if (res_occs.ket_occ(*site) != der_occs.ket_occ(*site))
                        compatible = false;
                    if (std::count(res_corr.site_vec().begin(), res_corr.site_vec().end(), *site) > 0)
                      bra_occs.set_occ(*site, res_occs.ket_occ(*site));
                    else if (std::count(der_corr.site_vec().begin(), der_corr.site_vec().end(), *site) > 0)
                      bra_occs.set_occ(*site, der_occs.ket_occ(*site));
                    else
                      throw pcps::Exception("expected an occupation for the cluster site");
                  }
                  if (compatible) {
                    coeff = -pcps::unity<S>() / der_corr[der_occs.cmpd()];
                    pcps::expectation_value<S>(slaves, slaves_used, jobs_in_progress, output_array.get(), cmpd_jac_index,
                                                    coeff, userinp, wfn, ham.ops().at(i), perturbation_level,
                                                    bra_occs, pcps::Cluster(), excluded);
                  }
                  if (myrank == 0 && userinp.verbose_print()) print_progress(total_ev, ++finished_ev, recorded_ev, sw);

                } // end loop over operators

              } // end loop over derivative occupations

            } // end loop over derivative correlators

          } // end loop over residual occupations

        } // end loop over residual correlators

        if (myrank == 0 && userinp.verbose_print()) {
          //std::cout << boost::format("%4i%%") % 100 << std::endl;
          std::cout << std::endl;
        }

      }

  };

  // create the worker that will compute the residual and jacobian precursers
  Local_Worker worker;

  // compute the residual and jacobian precursers in parallel
  const int n = wfn.nvars();
  boost::shared_array<S> data_array = pcps::compute_in_parallel(userinp, ham, wfn, perturbation_level, &worker, n*(n+3));

  // compute the energy
  const S energy = pcps::energy(userinp, ham, wfn, perturbation_level);

  // get MPI info
  const MPI::Comm & comm = MPI::COMM_WORLD;
  const int nproc = comm.Get_size();
  const int myrank = comm.Get_rank();

  // compute and apply the update step on the root process
  if (myrank == 0) {

    // The data array's size is nvars*(nvars+3).
    // The first nvars*nvars elements are the jacobian precursers
    // The next nvars elements are the energy derivative precursers
    // The next nvars elements are the residual precursers
    // The last nvars elements are the overlap matrix elements

    // Compute the residual using the energy, residual precursers, and overlap matrix.
    // This overwrites the residual precursers with the residual.
    {
      S * const overlap_elements = &data_array[n*n+2*n];
      S * const resid_precursers = &data_array[n*n+n];
      for (int i = 0; i < n; i++)
        resid_precursers[i] -= energy * overlap_elements[i];
    }

    // Compute the jacobian using the jacobian precursers, the energy derivative precursers, and the overlap elements.
    // This overwrites the jacobian precursers with the jacobian.
    {
      S * const jacob_precursers = &data_array[0];
      S * const energy_ders      = &data_array[n*n];
      S * const overlap_elements = &data_array[n*n+2*n];
      for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
          jacob_precursers[i*n+j] -= overlap_elements[i] * energy_ders[j];
    }

    // get pointers to the residual and jacobian
    S * const jacobian = &data_array[0];
    S * const residual = &data_array[n*n+n];

    // print the residual
    S residual_norm = 0.0;
    for (int i = 0; i < n; i++)
      residual_norm += pcps::conj(residual[i]) * residual[i];
    residual_norm = std::sqrt(residual_norm);
    if (userinp.verbose_print()) {
      std::cout << boost::format("residual norm = %12.4e") % std::abs(residual_norm) << std::endl;
      std::cout << std::endl;
      std::cout << "residual:" << std::endl;
      for (int i = 0; i < n; i++)
        std::cout << boost::format("%6i  %s") % i % pcps::format_number("%18.6f", residual[i]) << std::endl;
      std::cout << std::endl;
    }

    //// save the jacobian for error checking
    //boost::shared_array<S> saved_jac = pcps::allocate_shared_array<S>(n*n);
    //for (int i = 0; i < n*n; i++)
    //  saved_jac[i] = jacobian[i];

    // compute and print the singular values of the jacobian
    const int lwork = 5*n;
    boost::shared_array<S> u     = pcps::allocate_shared_array<S>(n*n);
    boost::shared_array<S> vt    = pcps::allocate_shared_array<S>(n*n);
    boost::shared_array<S> work  = pcps::allocate_shared_array<S>(lwork);
    boost::shared_array<double> s     = pcps::allocate_shared_array<double>(n);
    boost::shared_array<double> rwork = pcps::allocate_shared_array<double>(5*n);
    {
      int info;
      pcps::xgesvd('A', 'A', n, n, &jacobian[0], n, s.get(), vt.get(), n, u.get(), n, work.get(), lwork, rwork.get(), info);
      if (info != 0)
        throw pcps::Exception( (boost::format("xgesvd failed for jacobian with error code %i") % info).str() );
      //// check that the jacobian was computed correctly (jacobian = u * s * vt)
      //boost::shared_array<S> new_jac = pcps::allocate_shared_array<S>(n*n, pcps::zero<S>());
      //for (int i = 0; i < n; i++)
      //for (int j = 0; j < n; j++)
      //for (int k = 0; k < n; k++)
      //  new_jac[i*n+j] += u[i*n+k] * s[k] * vt[k*n+j];
      //for (int i = 0; i < n*n; i++)
      //  new_jac[i] -= saved_jac[i];
      //double max_error = 0.0;
      //for (int i = 0; i < n*n; i++)
      //  if (std::abs(new_jac[i]) > max_error)
      //    max_error = std::abs(new_jac[i]);
      //if (userinp.verbose_print())
      //  std::cout << boost::format("max error in new_jac = %.4e") % max_error << std::endl;
    }
    if (userinp.verbose_print()) {
      std::cout << std::endl;
      std::cout << "jacobian singular values:" << std::endl;
      for (int i = 0; i < n; i++)
        std::cout << boost::format("%6i  %s") % i % pcps::format_number("%18.6e", s[i]) << std::endl;
    }

    // compute the update vector
    const double shift = 0.0;
    boost::shared_array<S> update   = pcps::allocate_shared_array<S>(n, pcps::zero<S>());
    //for (int i = 0; i < n; i++)
    //for (int j = 0; j < n; j++)
    //for (int k = 0; k < n; k++)
    //  if (std::fabs(s[j]) > 1.0e-10)
    //    update[i] += vt[j*n+i] * (1.0 / ( s[j] + shift ) ) * u[k*n+j] * ( -1.0 * residual[k] );
    int m = 0;
    for (int i = 0; i < n && s[i] > 1.0e-10; i++, m++) {}
    boost::shared_array<S> intermed = pcps::allocate_shared_array<S>(m, pcps::zero<S>());
    pcps::xgemm('N', 'N', m, 1, n, pcps::unity<S>(), u.get(), n, residual, n, pcps::zero<S>(), intermed.get(), m);
    for (int i = 0; i < m; i++)
      intermed[i] = -intermed[i] / ( s[i] + shift );
    pcps::xgemm('N', 'N', n, 1, m, pcps::unity<S>(), vt.get(), n, intermed.get(), m, pcps::zero<S>(), update.get(), n);

//    // shift the jacobian's diagonal
//    const S shift = 0.01 * pcps::unity<S>();
//    for (int i = 0; i < n; i++)
//      jacobian[i*n + i] += shift;
//
//    // solve the linear equation ( J x = - r ) for the update vector
//    std::vector<S> update(residual.size(), pcps::zero<S>());
//    int info = 0;
//    std::vector<int> ipiv(n);
//    pcps::xgesv(n, 1, &jacobian[0], n, &ipiv[0], &residual[0], n, info);
//    if (info != 0)
//      throw pcps::Exception( (boost::format("xgesv failed for jacobian with error code %i") % info).str() );
//
//    // the negative of the update vector is now stored in residual

      // add the update step to the correlators

      // initialize an index for use in computing the compound residual index
      int base_index = 0;

      // loop over the correlators we are optimizing
      if (userinp.verbose_print())
        std::cout << std::endl;
      for (int c_ind = 0; c_ind < wfn.correlators().size(); c_ind++) {

        // get a reference to the correlator
        pcps::Correlator<S> & corr = wfn.correlators().at(c_ind);

        // fixed correlators and images of translationally invariant correlators are not updated
        if ( corr.fixed() || corr.ti() ) continue;

        // loop over the correlator's occupations
        for (pcps::Cluster bra_occs(userinp.nsites(), userinp.configs_per_site(), corr.site_vec()); !bra_occs.finished(); bra_occs++) {

          // get the residual element's compound index
          const int cmpd_res_index = wfn.corr_offsets()[c_ind] + bra_occs.cmpd();

          // update the correlator element
          corr[bra_occs.cmpd()] += ( 1.0 - damping ) * update[cmpd_res_index];

        } // end loop over occupations

        // normalize the elements
        corr.normalize();

        // print the new correlator
        if (userinp.verbose_print()) {
          pcps::print_correlator(corr);
          std::cout << std::endl;
        }

      } // end loop over correlators

  }

  // broadcast the new correlator elements to all processes
  comm.Barrier();
  for (typename std::vector< pcps::Correlator<S> >::iterator corr = wfn.correlators().begin(); corr != wfn.correlators().end(); corr++)
    comm.Bcast((void *)corr->data_ptr(), sizeof(S) * corr->nelements(), MPI::CHAR, 0);

}

// explicitly instantiate the template
template void pcps::projected_nr_update(const pcps::Input &,
                                        const pcps::Hamiltonian<double> &,
                                        pcps::Wavefunction<double> &,
                                        const int,
                                        const double);
template void pcps::projected_nr_update(const pcps::Input &,
                                        const pcps::Hamiltonian< std::complex<double> > &,
                                        pcps::Wavefunction< std::complex<double> > &,
                                        const int,
                                        const double);

// previous version without translational invariance or complex number support

////-------------------------------------------------------------------------------
//// pcps::projected_nr_update -- updates the correlators via a Newton-Raphson step
////-------------------------------------------------------------------------------
//
//template <class S>
//void pcps::projected_nr_update(const pcps::Input & userinp,
//                               const pcps::Hamiltonian<S> & ham,
//                               pcps::Wavefunction<S> & wfn,
//                               const int perturbation_level,
//                               const double damping) {
//
//  // not yet implemented for complex numbers
//  if (typeid(S) != typeid(double))
//    throw pcps::Exception("pcps::projected_nr_update is only implemented for real numbers.");
//
//  // define the worker class that will compute the residual and jacobian precursers in parallel
//  class Local_Worker : public pcps::Parallel_Worker_Base<S> {
//
//    public:
//
//      void do_work(const pcps::Input & userinp,
//                   const pcps::Hamiltonian<S> & ham,
//                   const pcps::Wavefunction<S> & wfn,
//                   const int perturbation_level,
//                   const boost::shared_array<S> & output_array,
//                   const int array_size,
//                   const int my_min,
//                   const int my_max,
//                   const std::vector<int> & slaves,
//                   int & slaves_used,
//                   int & jobs_in_progress)
//      {
//
//        // get MPI info
//        const MPI::Comm & comm = MPI::COMM_WORLD;
//        const int nproc = comm.Get_size();
//        const int myrank = comm.Get_rank();
//
//        // determine the number of variables in the wavefunction
//        const int nvars = wfn.nvars();
//
//        // The output array's size is nvars*(nvars+3).
//        // The first nvars*nvars elements are the jacobian precursers
//        // The next nvars elements are the energy derivative precursers
//        // The next nvars elements are the residual precursers
//        // The last nvars elements are the overlap matrix elements
//
//        // prepare a copy of the identity operator
//        const typename pcps::OpBase<S>::const_ptr identity_op(new pcps::OpIdentity<S>(pcps::unity<S>()));
//
//        // prepare an empty set of correlators to indicate that none should be excluded during expectation value calcualtions
//        typename pcps::Correlator<S>::Set excluded;
//
//        pcps::Stopwatch sw; sw.start();
//        double completion_counter = 0.0;
//
//        // initialize an index for use in computing the compound residual index
//        int base_res_index = 0;
//
//        // loop over the correlators that form the residual
//        for (int r_ind = 0; r_ind < wfn.correlators().size(); r_ind++) {
//
//          // get a reference to the residual correlator
//          const pcps::Correlator<S> & res_corr = wfn.correlators().at(r_ind);
//
//          // fixed correlators are not included
//          if (res_corr.fixed()) continue;
//
//          // loop over the correlator's occupations
//          for (pcps::Cluster res_occs(userinp.nsites(), userinp.configs_per_site(), res_corr.site_vec()); !res_occs.finished(); res_occs++) {
//
//            // get the residual element's compound index
//            const int cmpd_res_index = base_res_index + res_occs.cmpd();
//
//            const double old_completion = completion_counter;
//            while ( double(cmpd_res_index) / nvars > completion_counter ) { completion_counter += 0.01; }
//            if (myrank == 0 && userinp.verbose_print() && std::fabs(completion_counter - old_completion) > 1.0e-9) {
//              std::cout << boost::format("%4i%%  %10.4f elapsed seconds") % int(100 * completion_counter) % sw.elapsed_seconds() << std::endl;
//            }
//
//            // compute the overlap matrix element
//            if (myrank == 0) {
//              const int overlap_pos = nvars*nvars + 2*nvars + cmpd_res_index;
//              pcps::expectation_value<S>(slaves, slaves_used, jobs_in_progress, output_array.get(), overlap_pos, pcps::unity<S>(),
//                                              userinp, wfn, identity_op, perturbation_level, res_occs, pcps::Cluster(), excluded);
//            }
//
//            // compute the residual precurser element
//            const int res_precurser_pos = nvars*nvars + nvars + cmpd_res_index;
//            for (int i = my_min; i < my_max; i++)
//              pcps::expectation_value<S>(slaves, slaves_used, jobs_in_progress, output_array.get(), res_precurser_pos, pcps::unity<S>(),
//                                              userinp, wfn, ham.ops().at(i), perturbation_level, res_occs, pcps::Cluster(), excluded);
//
//            // compute the numerator derivative contribution to the energy derivative element (bra_occs = none, ket_occs = res_occs)
//            const int energy_der_pos = nvars*nvars + cmpd_res_index;
//            S coeff = pcps::unity<S>() / res_corr[res_occs.cmpd()];
//            for (int i = my_min; i < my_max; i++)
//              pcps::expectation_value<S>(slaves, slaves_used, jobs_in_progress, output_array.get(), energy_der_pos, coeff,
//                                              userinp, wfn, ham.ops().at(i), perturbation_level, pcps::Cluster(), res_occs, excluded);
//
//            // the denominator derivative contribution to the energy can be computed from
//            // the residual precursers and so does not need to be dealt with here
//
//            // initialize an index for use in computing the compound derivative index
//            int base_der_index = 0;
//
//            // loop over the correlators that we differentiate with respect to
//            for (int d_ind = 0; d_ind < wfn.correlators().size(); d_ind++) {
//
//              // get a reference to the derivative correlator
//              const pcps::Correlator<S> & der_corr = wfn.correlators().at(d_ind);
//
//              // fixed correlators are not included
//              if (der_corr.fixed()) continue;
//
//              // loop over the derivative correlator's occupations
//              for (pcps::Cluster der_occs(userinp.nsites(), userinp.configs_per_site(), der_corr.site_vec()); !der_occs.finished(); der_occs++) {
//                
//                // get the derivative element's compound index
//                const int cmpd_der_index = base_der_index + der_occs.cmpd();
//
//                // get the jacobian element's compound index
//                const int cmpd_jac_index = cmpd_res_index * wfn.nvars() + cmpd_der_index;
//
//                // compute the contribution from each hamiltonian operator
//                for (int i = my_min; i < my_max; i++) {
//
//                  // check if this operator contributes
//                  bool contributes = false;
//                  for (int j = 0; j < ham.ops().at(i)->sites().size(); j++)
//                    if (    ham.ops().at(i)->offsets().at(j) != 0
//                         && std::count(der_corr.site_vec().begin(), der_corr.site_vec().end(), ham.ops().at(i)->sites().at(j)) > 0)
//                      contributes = true;
//
//                  // if the operator does not contribute, skip it
//                  if (!contributes) continue;
//
//                  // compute the numerator contribution (bra_occs = res_occs, ket_occs = der_occs)
//                  coeff = pcps::unity<S>() / der_corr[der_occs.cmpd()];
//                  pcps::expectation_value<S>(slaves, slaves_used, jobs_in_progress, output_array.get(), cmpd_jac_index, coeff,
//                                                  userinp, wfn, ham.ops().at(i), perturbation_level, res_occs, der_occs, excluded);
//
//                  // compute the denominator contribution (bra_occs = res_occs + der_occs, ket_occs = none)
//                  std::set<int> bra_occ_sites;
//                  bra_occ_sites.insert(res_occs.sites().begin(), res_occs.sites().end());
//                  bra_occ_sites.insert(der_occs.sites().begin(), der_occs.sites().end());
//                  pcps::Cluster bra_occs(userinp.nsites(), userinp.configs_per_site(), bra_occ_sites);
//                  bool compatible = true;
//                  for (std::set<int>::const_iterator site = bra_occ_sites.begin(); site != bra_occ_sites.end() && compatible; site++) {
//                    if (    std::count(res_corr.site_vec().begin(), res_corr.site_vec().end(), *site) > 0
//                         && std::count(der_corr.site_vec().begin(), der_corr.site_vec().end(), *site) > 0 )
//                      if (res_occs.ket_occ(*site) != der_occs.ket_occ(*site))
//                        compatible = false;
//                    if (std::count(res_corr.site_vec().begin(), res_corr.site_vec().end(), *site) > 0)
//                      bra_occs.set_occ(*site, res_occs.ket_occ(*site));
//                    else if (std::count(der_corr.site_vec().begin(), der_corr.site_vec().end(), *site) > 0)
//                      bra_occs.set_occ(*site, der_occs.ket_occ(*site));
//                    else
//                      throw pcps::Exception("expected an occupation for the cluster site");
//                  }
//                  if (compatible) {
//                    pcps::Cluster ket_occs;
//                    coeff = -pcps::unity<S>() / der_corr[der_occs.cmpd()];
//                    pcps::expectation_value<S>(slaves, slaves_used, jobs_in_progress, output_array.get(), cmpd_jac_index, coeff,
//                                                    userinp, wfn, ham.ops().at(i), perturbation_level, bra_occs, ket_occs, excluded);
//                  }
//
//                } // end loop over operators
//
//              } // end loop over derivative occupations
//
//              // increase the base index by the number of variables in the correlator
//              base_der_index += der_corr.nelements();
//
//            } // end loop over derivative correlators
//
//          } // end loop over residual occupations
//
//          // increase the base index by the number of variables in the correlator
//          base_res_index += res_corr.nelements();
//
//        } // end loop over residual correlators
//
//        if (myrank == 0 && userinp.verbose_print()) {
//          std::cout << boost::format("%4i%%") % 100 << std::endl;
//          std::cout << std::endl;
//        }
//
//      }
//
//  };
//
//  // create the worker that will compute the residual and jacobian precursers
//  Local_Worker worker;
//
//  // compute the residual and jacobian precursers in parallel
//  const int n = wfn.nvars();
//  boost::shared_array<S> data_array = pcps::compute_in_parallel(userinp, ham, wfn, perturbation_level, &worker, n*(n+3));
//
//  // compute the energy
//  const S energy = pcps::energy(userinp, ham, wfn, perturbation_level);
//
//  // get MPI info
//  const MPI::Comm & comm = MPI::COMM_WORLD;
//  const int nproc = comm.Get_size();
//  const int myrank = comm.Get_rank();
//
//  // compute and apply the update step on the root process
//  if (myrank == 0) {
//
//    // The data array's size is nvars*(nvars+3).
//    // The first nvars*nvars elements are the jacobian precursers
//    // The next nvars elements are the energy derivative precursers
//    // The next nvars elements are the residual precursers
//    // The last nvars elements are the overlap matrix elements
//
//    // Compute the denominator energy derivative contribution using the residual precursers.
//    {
//      S * const resid_precursers = &data_array[n*n+n];
//      S * const energy_ders      = &data_array[n*n];
//      int base_res_index = 0;
//      for (int c_ind = 0; c_ind < wfn.correlators().size(); c_ind++) {
//        const pcps::Correlator<S> & corr = wfn.correlators().at(c_ind);
//        if (corr.fixed()) continue;
//        for (pcps::Cluster occs(userinp.nsites(), userinp.configs_per_site(), corr.site_vec()); !occs.finished(); occs++) {
//          const int res_index = base_res_index + occs.cmpd();
//          energy_ders[res_index] -= resid_precursers[res_index] / corr[occs.cmpd()];
//        }
//        base_res_index += corr.nelements();
//      }
//    }
//
//    // Compute the residual using the energy, residual precursers, and overlap matrix.
//    // This overwrites the residual precursers with the residual.
//    {
//      S * const overlap_elements = &data_array[n*n+2*n];
//      S * const resid_precursers = &data_array[n*n+n];
//      for (int i = 0; i < n; i++)
//        resid_precursers[i] -= energy * overlap_elements[i];
//    }
//
//    // Compute the jacobian using the jacobian precursers, the energy derivative precursers, and the overlap elements.
//    // This overwrites the jacobian precursers with the jacobian.
//    {
//      S * const jacob_precursers = &data_array[0];
//      S * const energy_ders      = &data_array[n*n];
//      S * const overlap_elements = &data_array[n*n+2*n];
//      for (int i = 0; i < n; i++)
//        for (int j = 0; j < n; j++)
//          jacob_precursers[i*n+j] -= overlap_elements[i] * energy_ders[j];
//    }
//
//    // get pointers to the residual and jacobian
//    S * const jacobian = &data_array[0];
//    S * const residual = &data_array[n*n+n];
//
//    // print the residual
//    S residual_norm = 0.0;
//    for (int i = 0; i < n; i++)
//      residual_norm += pcps::conj(residual[i]) * residual[i];
//    residual_norm = std::sqrt(residual_norm);
//    if (userinp.verbose_print()) {
//      std::cout << boost::format("residual norm = %12.4e") % std::abs(residual_norm) << std::endl;
//      std::cout << std::endl;
//      std::cout << "residual:" << std::endl;
//      for (int i = 0; i < n; i++)
//        std::cout << boost::format("%6i  %s") % i % pcps::format_number("%18.6f", residual[i]) << std::endl;
//      std::cout << std::endl;
//    }
//
//    //// save the jacobian for error checking
//    //boost::shared_array<S> saved_jac = pcps::allocate_shared_array<S>(n*n);
//    //for (int i = 0; i < n*n; i++)
//    //  saved_jac[i] = jacobian[i];
//
//    // compute and print the singular values of the jacobian
//    const int lwork = 5*n;
//    boost::shared_array<S> u     = pcps::allocate_shared_array<S>(n*n);
//    boost::shared_array<S> vt    = pcps::allocate_shared_array<S>(n*n);
//    boost::shared_array<S> work  = pcps::allocate_shared_array<S>(lwork);
//    boost::shared_array<double> s     = pcps::allocate_shared_array<double>(n);
//    boost::shared_array<double> rwork = pcps::allocate_shared_array<double>(5*n);
//    {
//      int info;
//      pcps::xgesvd('A', 'A', n, n, &jacobian[0], n, s.get(), vt.get(), n, u.get(), n, work.get(), lwork, rwork.get(), info);
//      if (info != 0)
//        throw pcps::Exception( (boost::format("xgesvd failed for jacobian with error code %i") % info).str() );
//      //// check that the jacobian was computed correctly (jacobian = u * s * vt)
//      //boost::shared_array<S> new_jac = pcps::allocate_shared_array<S>(n*n, pcps::zero<S>());
//      //for (int i = 0; i < n; i++)
//      //for (int j = 0; j < n; j++)
//      //for (int k = 0; k < n; k++)
//      //  new_jac[i*n+j] += u[i*n+k] * s[k] * vt[k*n+j];
//      //for (int i = 0; i < n*n; i++)
//      //  new_jac[i] -= saved_jac[i];
//      //double max_error = 0.0;
//      //for (int i = 0; i < n*n; i++)
//      //  if (std::abs(new_jac[i]) > max_error)
//      //    max_error = std::abs(new_jac[i]);
//      //if (userinp.verbose_print())
//      //  std::cout << boost::format("max error in new_jac = %.4e") % max_error << std::endl;
//    }
//    if (userinp.verbose_print()) {
//      std::cout << std::endl;
//      std::cout << "jacobian singular values:" << std::endl;
//      for (int i = 0; i < n; i++)
//        std::cout << boost::format("%6i  %s") % i % pcps::format_number("%18.6e", s[i]) << std::endl;
//    }
//
//    // compute the update vector
//    const double shift = 0.0;
//    boost::shared_array<S> update   = pcps::allocate_shared_array<S>(n, pcps::zero<S>());
//    //for (int i = 0; i < n; i++)
//    //for (int j = 0; j < n; j++)
//    //for (int k = 0; k < n; k++)
//    //  if (std::fabs(s[j]) > 1.0e-10)
//    //    update[i] += vt[j*n+i] * (1.0 / ( s[j] + shift ) ) * u[k*n+j] * ( -1.0 * residual[k] );
//    int m = 0;
//    for (int i = 0; i < n && s[i] > 1.0e-10; i++, m++) {}
//    boost::shared_array<S> intermed = pcps::allocate_shared_array<S>(m, pcps::zero<S>());
//    pcps::xgemm('N', 'N', m, 1, n, pcps::unity<S>(), u.get(), n, residual, n, pcps::zero<S>(), intermed.get(), m);
//    for (int i = 0; i < m; i++)
//      intermed[i] = -intermed[i] / ( s[i] + shift );
//    pcps::xgemm('N', 'N', n, 1, m, pcps::unity<S>(), vt.get(), n, intermed.get(), m, pcps::zero<S>(), update.get(), n);
//
////    // shift the jacobian's diagonal
////    const S shift = 0.01 * pcps::unity<S>();
////    for (int i = 0; i < n; i++)
////      jacobian[i*n + i] += shift;
////
////    // solve the linear equation ( J x = - r ) for the update vector
////    std::vector<S> update(residual.size(), pcps::zero<S>());
////    int info = 0;
////    std::vector<int> ipiv(n);
////    pcps::xgesv(n, 1, &jacobian[0], n, &ipiv[0], &residual[0], n, info);
////    if (info != 0)
////      throw pcps::Exception( (boost::format("xgesv failed for jacobian with error code %i") % info).str() );
////
////    // the negative of the update vector is now stored in residual
//
//      // add the update step to the correlators
//
//      // initialize an index for use in computing the compound residual index
//      int base_index = 0;
//
//      // loop over the correlators we are optimizing
//      if (userinp.verbose_print())
//        std::cout << std::endl;
//      for (int c_ind = 0; c_ind < wfn.correlators().size(); c_ind++) {
//
//        // get a reference to the correlator
//        pcps::Correlator<S> & corr = wfn.correlators().at(c_ind);
//
//        // fixed correlators are not being optimized
//        if (corr.fixed()) continue;
//
//        // loop over the correlator's occupations
//        for (pcps::Cluster bra_occs(userinp.nsites(), userinp.configs_per_site(), corr.site_vec()); !bra_occs.finished(); bra_occs++) {
//
//          // get the residual element's compound index
//          const int cmpd_res_index = base_index + bra_occs.cmpd();
//
//          // update the correlator element
//          corr[bra_occs.cmpd()] += ( 1.0 - damping ) * update[cmpd_res_index];
//
//        } // end loop over occupations
//
//        // normalize the elements
//        S max_value = pcps::zero<S>();
//        for (int i = 0; i < corr.nelements(); i++)
//          if (std::abs(corr[i]) > std::abs(max_value))
//            max_value = corr[i];
//        for (int i = 0; i < corr.nelements(); i++)
//            corr[i] = corr[i] / max_value;
//
//        // print the new correlator
//        if (userinp.verbose_print()) {
//          pcps::print_correlator(corr);
//          std::cout << std::endl;
//        }
//
//        // increase the base index by the number of variables in the correlator
//        base_index += corr.nelements();
//
//      } // end loop over correlators
//
//  }
//
//  // broadcast the new correlator elements to all processes
//  comm.Barrier();
//  for (typename std::vector< pcps::Correlator<S> >::iterator corr = wfn.correlators().begin(); corr != wfn.correlators().end(); corr++)
//    comm.Bcast((void *)corr->data_ptr(), sizeof(S) * corr->nelements(), MPI::CHAR, 0);
//
//}
//
//// explicitly instantiate the template
//template void pcps::projected_nr_update(const pcps::Input &,
//                                        const pcps::Hamiltonian<double> &,
//                                        pcps::Wavefunction<double> &,
//                                        const int,
//                                        const double);
//template void pcps::projected_nr_update(const pcps::Input &,
//                                        const pcps::Hamiltonian< std::complex<double> > &,
//                                        pcps::Wavefunction< std::complex<double> > &,
//                                        const int,
//                                        const double);
